% Estimation of the unshrunken model

%% Housekeeping
clc; clear all; close all;
workpath00 = pwd;

addpath('toolbox_xt');
addpath('toolbox_plot');

%% load data
load('data_all.mat');

%% Two panel variables
% xtflag and xtval: combofips
% xtflag2 and xtval2: year
% xtflag3 and xtval3: combofips/year
xtflag3 = xtflag + xtflag2/10000;
xtval3 = unique(xtflag3);

nobs = size(YY,1);
ncity = length(xtval);
nyear = 6;

%% Data for the esitmation
% sel_est = true(nobs,1);  %include all

[~,~,info] = xt_reg(YY, ones(nobs,1), xtflag3, xtval3);
sel_est = (info.n>=5);


est.YY  = YY(sel_est,:);
est.xtflag = xtflag(sel_est,:);
est.xtval  = unique(est.xtflag);
est.xtflag2 = xtflag2(sel_est,:);
est.xtval2  = unique(est.xtflag2);
est.xtflag3 = xtflag3(sel_est,:);
est.xtval3  = unique(est.xtflag3);
est.nobs    = size(xtflag,1);

% Other controls
YD = [year_dummy2, year_dummy3, year_dummy4, year_dummy5, year_dummy6];

% estimation
est.YD = YD(sel_est,:);
est.XX = XX(sel_est,:);
est.logDC = logDC(sel_est,:);
est.DD = DD(sel_est,:);

est.ZZ = ZZ(sel_est,:);
[~,temp_ZZ2] = xt_reg(est.ZZ(:,2), ones(size(est.ZZ,1),1), est.xtflag, est.xtval);
est.ZZ2 = [ones(size(temp_ZZ2,1),1), temp_ZZ2];

estinfo.islog = 2;
estinfo.FitMethod = 'REML';

% unshrunken
% est_us_m2 = estimator_us_p1(est, estinfo);
% est_us_m2.model.MSE

est_us_m5 = estimator_us_p4s3(est, estinfo); %almost impossible to do this esitmation...
est_us_m5.model.MSE

% shrinkage
% est_s_m2 = estimator_s_A_p1s3(est, estinfo);
est_s_m5 = estimator_s_A_p4s3(est, estinfo); %

% save
save('results_shrunken_s3.mat');


